This report runs base, AR1, and AR2 models (variable ~ time) for annual SEV meteorological data using the nlme R package. Mean annual air temperature, minimum annual air temperature, maximum annual air temperature, and total annual precipitation are modeled. More sophisticated models will likely be required after this preliminary analysis.
library(tidyverse)
library(lubridate)
library(nlme)
library(MuMIn)
library(gridExtra)
library(kableExtra)
# Load yearly data ---------------------------------------------
path_to_files <- "../data/processed_data/"
# load and only keep vars of interest
y <- read_csv(paste0(path_to_files, "met_yearly_gap_filled.csv")) %>%
mutate(sta = as.factor(sta)) %>%
select(sta:ppt)
str(y)
## tibble [110 × 6] (S3: tbl_df/tbl/data.frame)
## $ sta : Factor w/ 4 levels "40","42","49",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ year : num [1:110] 1988 1989 1990 1991 1992 ...
## $ airt : num [1:110] 13.8 13.8 13.8 12.5 12.3 ...
## $ maxair: num [1:110] 35 35 35 37.7 38.6 ...
## $ minair: num [1:110] -7.94 -7.94 -7.94 -13.76 -15.11 ...
## $ ppt : num [1:110] 151 140 232 335 241 ...
summary(y)
## sta year airt maxair minair
## 40:34 Min. :1988 Min. :11.91 Min. :32.01 Min. :-29.930
## 42:33 1st Qu.:2001 1st Qu.:13.15 1st Qu.:37.13 1st Qu.:-16.495
## 49:23 Median :2008 Median :13.82 Median :38.70 Median :-14.780
## 50:20 Mean :2007 Mean :13.89 Mean :38.50 Mean :-14.900
## 3rd Qu.:2015 3rd Qu.:14.67 3rd Qu.:39.84 3rd Qu.:-12.982
## Max. :2021 Max. :16.43 Max. :48.72 Max. : -3.797
## ppt
## Min. :122.0
## 1st Qu.:220.9
## Median :247.9
## Mean :262.5
## 3rd Qu.:301.7
## Max. :599.4
Add a time variable for nlme modeling and split data into individual met stations.
# function to select station and variable and add 'time' var
prepare_data_for_nlme <- function(data, sta) {
# data = dataframe to use
# sta = met station (in quotes)
data %>%
filter(sta == {{ sta }}) %>%
mutate(time = as.numeric(as.factor(year)))
}
m40 <- prepare_data_for_nlme(y, "40")
m42 <- prepare_data_for_nlme(y, "42")
m49 <- prepare_data_for_nlme(y, "49")
m50 <- prepare_data_for_nlme(y, "50")
This code produces several functions that will be useful for producing cleaner code.
# function for plotting preliminary yearly graph -----------------------------------------------------------
plot_yearly_prelim <- function(data, var, title, y_axis_label) {
# data = dataset
# var = variable (not in quotes)
# title = title for graph (in quotes)
# y_axis_label = label for y-axis (in quotes)
ggplot(data, aes(x = year, y = {{ var }}, col = sta)) +
geom_line(color = "burlywood") +
geom_smooth(method = "lm", color = "black", size = 0.6) +
facet_wrap(~ sta) +
labs(title = title,
x = "Year",
y = y_axis_label)
}
# function for plotting graphs with best nlme model results ------------------------------------------------
plot_yearly_results <- function(data, var, coeffs, title, y_axis_label) {
# data = dataset
# var = variable (not in quotes)
# coeffs = object containing coefficients and p-values of best model
# title = title for graph (in quotes)
# y_axis_label = label for y-axis (in quotes)
ggplot(data, aes(x = year, y = {{ var }})) +
geom_line(color = "burlywood") +
geom_smooth(method = "lm", color = "black", size = 0.6) +
labs(title = title,
x = "Year",
y = y_axis_label,
caption = paste("slope = ", round(coeffs[2,1], 3), "\n(p-value = ", coeffs[2,2], ")"))
}
# functions to run nlme models ----------------------------------------------------------
# base model function
base_nlme_model <- function(data, var) {
# data = dataframe
# var = var to model
model_specification = as.formula(paste0(var, " ~ time"))
gls(model = model_specification,
data = data,
method = "ML",
na.action = na.omit)
}
# AR1 model function
AR1_nlme_model <- function(data, var) {
# data = dataframe
# var = var to model
model_specification = as.formula(paste0(var, " ~ time"))
gls(model = model_specification,
data = data,
correlation = corARMA(form = ~ time, p = 1),
method = "ML",
na.action = na.omit)
}
# AR2 model function
AR2_nlme_model <- function(data, var) {
# data = dataframe
# var = var to model
model_specification = as.formula(paste0(var, " ~ time"))
gls(model = model_specification,
data = data,
correlation = corARMA(form = ~ time, p = 2),
method = "ML",
na.action = na.omit)
}
plot_yearly_prelim(y, airt, "Annual Mean Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
# Met 40:
m40_mbase <- base_nlme_model(m40, "airt")
m40_mAR1 <- AR1_nlme_model(m40, "airt")
m40_mAR2 <- AR2_nlme_model(m40, "airt")
AICc scores:
AICc(m40_mbase, m40_mAR1, m40_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_mbase | 3 | 53.68867 |
| m40_mAR1 | 4 | 55.62003 |
| m40_mAR2 | 5 | 58.37968 |
summary(m40_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 52.88867 57.46775 -23.44433
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 13.271801 0.17431449 76.13711 0.0000
## time 0.021443 0.00868861 2.46794 0.0191
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.1560729 -0.5446061 0.2023188 0.6881964 1.6138251
##
## Residual standard error: 0.4821986
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, airt, m40_coeffs, "Met 40 - Deep Well \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 42:
m42_mbase <- base_nlme_model(m42, "airt")
m42_mAR1 <- AR1_nlme_model(m42, "airt")
m42_mAR2 <- AR2_nlme_model(m42, "airt")
AICc scores:
AICc(m42_mbase, m42_mAR1, m42_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_mbase | 3 | 50.38362 |
| m42_mAR1 | 4 | 52.94158 |
| m42_mAR2 | 5 | 55.65049 |
summary(m42_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 49.55604 54.04556 -21.77802
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 12.618607 0.17205515 73.34048 0.0000
## time 0.011535 0.00883012 1.30629 0.2011
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.987976067 -0.644320891 -0.002985058 0.916472339 1.576226874
##
## Residual standard error: 0.468135
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, airt, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 49:
m49_mbase <- base_nlme_model(m49, "airt")
m49_mAR1 <- AR1_nlme_model(m49, "airt")
m49_mAR2 <- AR2_nlme_model(m49, "airt")
AICc scores:
AICc(m49_mbase, m49_mAR1, m49_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_mbase | 3 | 34.60607 |
| m49_mAR1 | 4 | 37.16381 |
| m49_mAR2 | 5 | 39.81788 |
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 33.34291 36.74939 -13.67145
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 14.203751 0.19776763 71.82041 0.0000
## time 0.024676 0.01442369 1.71081 0.1018
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.1832629 -0.6706411 0.2109746 0.6989873 1.7965345
##
## Residual standard error: 0.4384421
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, airt, m49_coeffs, "Met 49 - Five Points \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 50:
m50_mbase <- base_nlme_model(m50, "airt")
m50_mAR1 <- AR1_nlme_model(m50, "airt")
m50_mAR2 <- AR2_nlme_model(m50, "airt")
AICc scores:
AICc(m50_mbase, m50_mAR1, m50_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_mbase | 3 | 33.37573 |
| m50_mAR1 | 4 | 35.87783 |
| m50_mAR2 | 5 | 39.37225 |
summary(m50_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 31.87573 34.86293 -12.93786
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 14.831383 0.22625535 65.55152 0.0000
## time 0.051514 0.01888743 2.72741 0.0138
##
## Correlation:
## (Intr)
## time -0.877
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.81172625 -0.64841217 -0.09216346 0.63402578 1.96854695
##
## Residual standard error: 0.462067
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, airt, m50_coeffs, "Met 50 - Blue Grama \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# airt results
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
plot_yearly_prelim(y, minair, "Annual Minimum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
# Met 40:
m40_mbase <- base_nlme_model(m40, "minair")
m40_mAR1 <- AR1_nlme_model(m40, "minair")
m40_mAR2 <- AR2_nlme_model(m40, "minair")
AICc scores:
AICc(m40_mbase, m40_mAR1, m40_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_mbase | 3 | 189.6589 |
| m40_mAR1 | 4 | 191.1119 |
| m40_mAR2 | 5 | 193.7564 |
summary(m40_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 188.8589 193.438 -91.42944
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -13.656159 1.2874555 -10.607092 0.0000
## time -0.136817 0.0641725 -2.132025 0.0408
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.6474716 -0.3773667 0.0234212 0.4512047 1.7215955
##
## Residual standard error: 3.561432
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, minair, m40_coeffs, "Met 40 - Deep Well \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 42:
m42_mbase <- base_nlme_model(m42, "minair")
m42_mAR1 <- AR1_nlme_model(m42, "minair")
m42_mAR2 <- AR2_nlme_model(m42, "minair")
AICc scores:
AICc(m42_mbase, m42_mAR1, m42_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_mbase | 3 | 189.4901 |
| m42_mAR1 | 4 | 190.6799 |
| m42_mAR2 | 5 | 192.9317 |
summary(m42_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 188.6625 193.152 -91.33124
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -13.220095 1.4158546 -9.337184 0.0000
## time -0.094695 0.0726637 -1.303196 0.2021
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.52824489 -0.32599081 -0.05008804 0.56078847 2.47065756
##
## Residual standard error: 3.852318
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, minair, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 49:
m49_mbase <- base_nlme_model(m49, "minair")
m49_mAR1 <- AR1_nlme_model(m49, "minair")
m49_mAR2 <- AR2_nlme_model(m49, "minair")
AICc scores:
AICc(m49_mbase, m49_mAR1, m49_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_mbase | 3 | 123.2896 |
| m49_mAR1 | 4 | 125.6878 |
| m49_mAR2 | 5 | 128.9351 |
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 122.0264 125.4329 -58.01321
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -14.590316 1.3596668 -10.730802 0.00
## time 0.006294 0.0991639 0.063475 0.95
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.0179895 -0.3267865 0.1783908 0.5709077 1.3371903
##
## Residual standard error: 3.014321
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, minair, m49_coeffs, "Met 49 - Five Points \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 50:
m50_mbase <- base_nlme_model(m50, "minair")
m50_mAR1 <- AR1_nlme_model(m50, "minair")
m50_mAR2 <- AR2_nlme_model(m50, "minair")
AICc scores:
AICc(m50_mbase, m50_mAR1, m50_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_mbase | 3 | 112.1543 |
| m50_mAR1 | 4 | 115.2297 |
| m50_mAR2 | 5 | 118.1863 |
summary(m50_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 110.6543 113.6415 -52.32717
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -14.264158 1.6215367 -8.796691 0.0000
## time 0.072444 0.1353633 0.535179 0.5991
##
## Correlation:
## (Intr)
## time -0.877
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.8049348 -0.3200412 0.1852600 0.5306493 1.1520424
##
## Residual standard error: 3.311562
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, minair, m50_coeffs, "Met 50 - Blue Grama \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# minair results
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
plot_yearly_prelim(y, maxair, "Annual Maximum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
# Met 40:
m40_mbase <- base_nlme_model(m40, "maxair")
m40_mAR1 <- AR1_nlme_model(m40, "maxair")
m40_mAR2 <- AR2_nlme_model(m40, "maxair")
AICc scores:
AICc(m40_mbase, m40_mAR1, m40_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_mbase | 3 | 140.3235 |
| m40_mAR1 | 4 | 135.6487 |
| m40_mAR2 | 5 | 136.8874 |
The AR1 model has the lowest AICc score. This model is used for plotting.
summary(m40_mAR1)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 134.2694 140.3748 -63.13468
##
## Correlation Structure: AR(1)
## Formula: ~time
## Parameter estimate(s):
## Phi
## 0.4530549
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 37.64031 0.9617813 39.13604 0.0000
## time 0.08334 0.0473790 1.75903 0.0881
##
## Correlation:
## (Intr)
## time -0.862
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6634576 -0.5746671 -0.1493423 0.6655313 3.1746966
##
## Residual standard error: 1.732308
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mAR1)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, maxair, m40_coeffs, "Met 40 - Deep Well \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 42:
m42_mbase <- base_nlme_model(m42, "maxair")
m42_mAR1 <- AR1_nlme_model(m42, "maxair")
m42_mAR2 <- AR2_nlme_model(m42, "maxair")
AICc scores:
AICc(m42_mbase, m42_mAR1, m42_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_mbase | 3 | 129.8644 |
| m42_mAR1 | 4 | 127.7573 |
| m42_mAR2 | 5 | 130.0821 |
The AR1 model has the lowest AICc score. This model is used for plotting.
summary(m42_mAR1)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 126.3287 132.3147 -59.16435
##
## Correlation Structure: AR(1)
## Formula: ~time
## Parameter estimate(s):
## Phi
## 0.3867261
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 34.94019 0.8290183 42.14647 0.0000
## time 0.06777 0.0421570 1.60763 0.1181
##
## Correlation:
## (Intr)
## time -0.864
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.94732099 -0.54279907 -0.08489296 0.66776673 2.67346131
##
## Residual standard error: 1.572183
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mAR1)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, maxair, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 49:
m49_mbase <- base_nlme_model(m49, "maxair")
m49_mAR1 <- AR1_nlme_model(m49, "maxair")
m49_mAR2 <- AR2_nlme_model(m49, "maxair")
AICc scores:
AICc(m49_mbase, m49_mAR1, m49_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_mbase | 3 | 116.1389 |
| m49_mAR1 | 4 | 117.1260 |
| m49_mAR2 | 5 | 116.5718 |
The base model has the lowest AICc and is used for plotting.
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 114.8758 118.2822 -54.43788
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 41.07945 1.1639166 35.29415 0.0000
## time -0.10151 0.0848873 -1.19584 0.2451
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.3311794 -0.5068364 -0.1698531 0.1037714 3.2757736
##
## Residual standard error: 2.580352
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, maxair, m49_coeffs, "Met 49 - Five Points \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 50:
m50_mbase <- base_nlme_model(m50, "maxair")
m50_mAR1 <- AR1_nlme_model(m50, "maxair")
m50_mAR2 <- AR2_nlme_model(m50, "maxair")
AICc scores:
AICc(m50_mbase, m50_mAR1, m50_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_mbase | 3 | 63.96627 |
| m50_mAR1 | 4 | 63.49305 |
| m50_mAR2 | 5 | 66.10862 |
The AR1 model has the lowest AICc score. This model is used for plotting.
summary(m50_mAR1)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 60.82638 64.80931 -26.41319
##
## Correlation Structure: AR(1)
## Formula: ~time
## Parameter estimate(s):
## Phi
## -0.4018488
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 38.90612 0.3271151 118.93708 0.0000
## time 0.07382 0.0274879 2.68568 0.0151
##
## Correlation:
## (Intr)
## time -0.882
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.65156196 -0.77367285 -0.05306683 0.49518631 2.02742794
##
## Residual standard error: 0.9854843
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mAR1)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, maxair, m50_coeffs, "Met 50 - Blue Grama \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# maxair results
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
plot_yearly_prelim(y, ppt, "Total Annual Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'
# Met 40:
m40_mbase <- base_nlme_model(m40, "ppt")
m40_mAR1 <- AR1_nlme_model(m40, "ppt")
m40_mAR2 <- AR2_nlme_model(m40, "ppt")
AICc scores:
AICc(m40_mbase, m40_mAR1, m40_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_mbase | 3 | 376.9800 |
| m40_mAR1 | 4 | 379.0143 |
| m40_mAR2 | 5 | 380.8387 |
summary(m40_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 376.18 380.7591 -185.09
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 239.73384 20.234522 11.847764 0.0000
## time -0.79358 1.008579 -0.786827 0.4372
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.7633224 -0.6843958 0.0941523 0.8117423 1.7513612
##
## Residual standard error: 55.97387
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, ppt, m40_coeffs, "Met 40 - Deep Well \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 42:
m42_mbase <- base_nlme_model(m42, "ppt")
m42_mAR1 <- AR1_nlme_model(m42, "ppt")
m42_mAR2 <- AR2_nlme_model(m42, "ppt")
AICc scores:
AICc(m42_mbase, m42_mAR1, m42_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_mbase | 3 | 386.6521 |
| m42_mAR1 | 4 | 389.0729 |
| m42_mAR2 | 5 | 390.7934 |
summary(m42_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 385.8245 390.314 -189.9122
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 391.6703 28.079405 13.948670 0.0000
## time -3.5405 1.441075 -2.456826 0.0198
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.60096602 -0.53817435 -0.06228751 0.56988895 3.13557611
##
## Residual standard error: 76.39964
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, ppt, m42_coeffs, "Met 42 - Cerro Montoso \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 49:
m49_mbase <- base_nlme_model(m49, "ppt")
m49_mAR1 <- AR1_nlme_model(m49, "ppt")
m49_mAR2 <- AR2_nlme_model(m49, "ppt")
AICc scores:
AICc(m49_mbase, m49_mAR1, m49_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_mbase | 3 | 253.8175 |
| m49_mAR1 | 4 | 256.5620 |
| m49_mAR2 | 5 | 259.8661 |
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 252.5543 255.9608 -123.2772
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 254.18442 23.215104 10.949096 0.0000
## time -1.99381 1.693136 -1.177584 0.2521
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.01258367 -0.67241549 0.07653863 0.62248810 2.09072501
##
## Residual standard error: 51.46686
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, ppt, m49_coeffs, "Met 49 - Five Points \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
# Met 50:
m50_mbase <- base_nlme_model(m50, "ppt")
m50_mAR1 <- AR1_nlme_model(m50, "ppt")
m50_mAR2 <- AR2_nlme_model(m50, "ppt")
AICc scores:
AICc(m50_mbase, m50_mAR1, m50_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_mbase | 3 | 222.6721 |
| m50_mAR1 | 4 | 225.6402 |
| m50_mAR2 | 5 | 229.2037 |
summary(m50_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 221.1721 224.1593 -107.586
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 245.26421 25.695541 9.545011 0.0000
## time 0.26436 2.145022 0.123244 0.9033
##
## Correlation:
## (Intr)
## time -0.877
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.95082391 -0.49761448 -0.02635629 0.88260113 1.46242374
##
## Residual standard error: 52.47638
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, ppt, m50_coeffs, "Met 50 - Blue Grama \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
# ppt results
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'